Event-driven simulations of a plastic, spiking neural network 
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We consider a fully-connected network of leaky integrate-and-fire neurons with spike-timing- 
dependent plasticity. The plasticity is controlled by a parameter representing the expected weight of 
a synapse between neurons that are firing randomly with the same mean frequency. For low values 
of the plasticity parameter, the activities of the system are dominated by noise, while large values 
of the plasticity parameter lead to self-sustaining activity in the network. We perform event-driven 
simulations on finite-size networks with up to 128 neurons to find the stationary synaptic weight 
conformations for different values of the plasticity parameter. In both the low and high activity 
regimes, the synaptic weights are narrowly distributed around the plasticity parameter value consis- 
tent with the predictions of mean-field theory. However, the distribution broadens in the transition 
region between the two regimes, representing emergent network structures. Using a pseudophysi- 
cal approach for visualization, we show that the emergent structures are of "path" or "hub" type, 
observed at different values of the plasticity parameter in the transition region. 

PACS numbers: 87.18.Sn, 87.19.lj, 87.19. Iw 



I. INTRODUCTION 

Neurons can form plastic networks through connecting 
synapses with weights that are changed dynamically by 
the neural activity in the form of neural spikes. It has 
been established that the precise timing of these spikes 
determines whether and how the synaptic weights will 
be increased (potentiated) or decreased (depressed) [1- 
4] . With better understanding of spike-timing-dependent 
plasticity (STDP), it becomes increasingly important 
to find out its implications on the underlying network 
structure and, consequently, on the neural activity itself. 
Theoretical models of a plastic neural network typically 
consist of three components: neural dynamics, synap- 
tic transmission, and network plasticity [5, G]. As has 
been shown previously [6] , a simple self-consistent mean- 
field scheme can be constructed when these three ingre- 
dients of a plastic neural network are given. Reducing 
the neural and synaptic dynamics to a determination of 
response functions characterizing a mean-field medium, 
the network plasticity is left governed by one-dimensional 
random-walk dynamics. Such a mean-field scheme was 
applied [G] to a network of leaky integrate-and-fire neu- 
rons [7-9], coupled through a jump-and-dccay synaptic 
conductance dynamics [10], and under STDP rules [2], 
controlled by a plasticity parameter w* representing the 
expected value of synaptic weights when firings of the 
neurons are purely driven by noise. While we do not 
expect that there is a broad range of cell types with 
different plasticity parameter values but otherwise sim- 
ilar physiological characteristics, different values of w* 
can correspond to different stages of a (quasi-statically) 
growing and developing network. The mean-field theory 
(MFT) predicts a first-order phase transition and hys- 
teresis from a low w* regime of noise-driven activity to 
a high w* regime of self-sustaining activity. It also pre- 
dicts a narrow synaptic weight distribution as long as 
the overall rate of change of synaptic strength is low. 



In the current study, we perform intensive event-driven 
simulations [11] on the same model network to compare 
with the predictions of MFT; the simulations also reveal 
emergent network dynamics and structures that are not 
captured by the MFT. 

In Section II that follows, we briefly describe the simu- 
lated model of a plastic neural network, which was elab- 
orated in our MFT study [()] . We then present in Section 
III the method and algorithm of the event-driven simu- 
lations used to explore the dynamics of the model. The 
results of the simulations for static uniform networks are 
reported in Section IV, where we present the mean ac- 
tivity levels for networks of different sizes as functions 
of synaptic input. These activity levels agree well with 
the MFT predictions in the stable phases of low and high 
activity but show a gradual increase in a transition re- 
gion (w* < w* < W2, from the onset to the persistence 
of threshold flring events) instead of the sharp first-order 
jump predicted by MFT. We also analyze the sporadic 
patterns of self-sustaining threshold firings in the transi- 
tion region to identify two types of threshold firing events. 
In Section V, we present simulation results with high res- 
olution in the plasticity parameter w* for a plastic net- 
work with = 32 neurons that has reached a station- 
ary state of plasticity. The change in the synaptic- weight 
conformation of the network in the transition region man- 
ifests itself in the synaptic weight distribution, which is 
seen to broaden twice, along with a bimodal elevation of 
the average firing activity compared to that of a static 
network. To characterize the emergent synaptic- weight 
structure of the network in the transition region, we em- 
ploy a pseudophysical approach for visualization in Sec- 
tion VI to generate a 2D layout of the network. Such an 
approach reveals a path or loop conformation near the 
lower end of the transition region {w* ^ wi) for network 
sizes up to « 32 and a hub conformation near the 
high end {w* ~ ^2) for N fa 24 and greater. We then 
conclude and summarize our findings in Section VII. 
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II. MODEL 

While our method is apphcable to other combinations 
of models of neurons, synapses and plasticity, we follow 
the same choices made in our mean- field approach [(>], 
which we describe briefly below. 

In the integrate-and-fire model for the neurons, the 
state of a neuron i is described by a membrane poten- 
tial Vi, which follows the differential equation of a leaky 
integrator [12] 



Table I. Values of parameters used in calculations 



integrate and fire 



resting potential Vb 

leak time 
firing threshold Vih 
reset potential Vi: 
reversal potential TZ: 



-55 mV 
20 ms 
-54 mV 
-80 mV 
mV 



TUM model 

decay time Tr>- 20 ms 
recovery time tr: 200 ms 
release fraction u: 0.5 
noise frequency Ajv: 1 Hz 
plasticity rate r: 0.01 



to the active state by each presynaptic spike. With the 
conservation rule 



' dt 



= Vo-v, + G,in-v,: 



(1) 



(5) 



where is the leak time for the membrane charge, Vq is 
the resting potential when the neuron is in the quiescent 
state, and TZ is the reversal potential for the ion channels 
on the synapses. The total synaptic conductance Gi for 
the neuron is given by the sum 



(2) 



over all presynaptic neurons, j. The synaptic weights 
Wj^i define the network and the same active transmitter 
fraction Yj is assumed for all efferent synapses of neuron 
j. In addition to the continuous dynamics (1), a neuron 
fires when its membrane potential reaches a threshold 
value, Vth- Then, its membrane potential drops imme- 
diately to a reset value, V,.. The action potential of the 
integrate-and-fire model is assumed to be instantaneous 
and is not modeled explicitly. The spike train produced 
by the neuron i is defined as the function 



(3) 



where ti^n is the time when the neuron i fires for the n-th 
time. 

The fraction Yj of the active transmitters is described 
by the Tsodyks-Uziel-Markram (TUM) model [10] of 
neural transmission, where the transmitters are dis- 
tributed in three states: "active", with the fraction Y; 
"inactive", with the fraction Z; and "ready-to-release", 
with the fraction X. For efferent synapses of a presynap- 
tic neuron j, these fractions follow the dynamics [10] 



dt — ' ■ 
dY^ 
dt 
dZj 
dt 



TR 

Y, 

+ uSjX-j 

td 



(4) 



TD 



TR 



where td is the decay time of active transmitters to the 
inactive state, tr is the recovery time for the inactive 
transmitters to the ready-to-release state, and u is the 
fraction of ready-to-release transmitters that is released 



there are two independent variables per presynaptic neu- 
ron. Consistent with the TUM dynamics, the values of 
the factors multiplying Sj at the discontinuities arc to be 
evaluated immediately before the discontinuities. 

While the integrate-and-fire and TUM dynamics are 
both deterministic, we model the stochasticity of the net- 
work with additional noise-driven firing events following 
Poisson statistics with the frequency Xjsi for each neu- 
ron. The noise-driven firings are treated the same way 
as threshold firings, that is, the membrane potentials are 
brought instantaneously to the reset value and the 
firing times are included in the spike trains (3) of the 
transmitter dynamics (4). 

To minimize the computational cost, the active trans- 
mitter fractions Yj double as the exponentially decaying 
"window functions" in our version of the plasticity rules 
suggested by van Rossum et al. [2] . The synaptic weights 
are taken to follow the dynamical equation 



(6) 



where A is the parameter for additive potentiation and r 
is the parameter for multiplicative depression. We define 
the plasticity parameter as the ratio 



, A 

w = —, 



(7) 



which is equal to the expectation value of the synaptic 
weight when we have the symmetry (YjSi) ~ {YiSj)^ e.g., 
when Si and Sj are both Poisson spike trains of the same 
frequency. Fixing u;* leaves the depression factor r as an 
overall control parameter for the rate of plasticity. 

It is common for a model of a biological system to 
carry a large number of empirical parameters. Instead of 
exploring all possible ranges of these parameters, we fix 
them with physiologically plausible values that are com- 
monly found in the literature. Unless otherwise stated, 
the values of the parameters are as listed in Table I. 



III. SIMULATION METHOD 

The disparity between the time scales of spiking activ- 
ities and neural plasticity posts a significant challenge to 
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computer simulations of plastic neural networks. Partic- 
ularly for STDP, where precise timing is crucial in deter- 
mining synaptic changes, we can not alleviate the com- 
putation requirement through use of larger integration 
time steps. However, for continuous dynamics that can 
be solved analytically, one can improve the efficiency of 
computation through an event- driven approach similar 
to that described by Brette [()] , which gives machine pre- 
cision timing for the spikes and requires limited amount 
of computation upon each spike production. 

In the event-driven approach, instead of calculating the 
state of the system for a fixed increment of time, one cal- 
culates the time for the next discontinuous event. This 
approach is feasible when the continuous dynamics of the 
system is solvable so that the time for the next discon- 
tinuous event can be evaluated efficiently. For the leaky 
integratc-and-fire model considered, an analytical form 
can be written down for the trajectory of the normalized 
membrane potential [!)] 



10"^ r 



V^-Va 

n-Vn 



El {x) - El (G) + 



GeG 



(8) 



where Ei is the exponential integral function of the first 
kind, X = Ge"*/"^'", and v = Vi (0), G = G; (0) are the 
initial states of the neuron. For the TUM dynamics, the 
transmitter fractions follow the trajectories 



Z, 



-t/TD 



Y 
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TR 
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TR 



TR - TD 



-Y^ (9) 



where Y = Yi (0) and Z = Zi (0) are their initial values. 
With the trajectory (8), one can solve for the time of the 
next threshold crossing Vi — > Vth through, e.g., a root 
solver. The computation time for a root solver to find 
the position of the crossing depends only logarithmically 
on the precision required, compared to the linear increase 
for the fixed-time-step approach. Further improvement 
of the computational efficiency can be achieved by, e.g., 
precomputing a lookup table for the time to threshold 
crossing given the current state of a neuron and interpo- 
lating when the desired value falls between prccomputcd 
points. For the simulated network with Poisson noise, 
the algorithm is as outlined below: 



1. For each neuron i in the system, calculate the time- 
to-fire tl^f from the continuous, deterministic dy- 
namics of the model. 



2. Draw the time to its next noise-triggered firing 
from an exponential distribution with the expecta- 
tion value for each neuron i. 

3. Find the minimum of the times tj^f and t]^ among 
all neurons to be the time to the next discrete event 
At of the system. 

4. Advance the simulated clock by At and update the 
state of the system using Eqs. (8) and (9). 
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Figure 1. (Color online) Mean activity level vs. total synap- 
tic input (Kw) for fully-connected uniform networks (node 
degree A" = A*' — 1) with static synaptic weight (w). Dashed 
lines are predictions from mean-field theory without fluctua- 
tion corrections (long dashes) and with the inclusion of shot- 
noise-like fluctuations (short dashes) in total synaptic con- 
ductance for a network with K — 63 afferent synapses per 
neuron. 



5. Fire the selected neuron i by reseting its membrane 
potential to Vi = Vr, and increasing the affected 
active transmitter fraction Yi by w (1 — 1^ — Zi). 

6. Repeat this process from 1. 

For the current model, the simulated dynamics boils 
down to the evaluation of the time-to-fire tttf {V, G) given 
the initial membrane potential V and total synaptic con- 
ductance G as well as the evaluation of the trajectories 
(8) and (9). We note that instead of the common ap- 
proach of defining the model with the set of differential 
equations (1) and (4), it is equally valid and computa- 
tionally preferable to define the continuous dynamics of 
the model with the trajectories (8) and (9) for the neu- 
rons and synapses. 



IV. STATIC NETWORK 

Before incorporating plasticity, we first perform event- 
driven simulations on fully-connected networks of up to 
128 neurons with uniform, fixed synaptic weights w. This 
exercise will prepare us for the inclusion of plasticity and 
also reveal the emergence of "structure" in the network. 
The average activity levels for different system sizes are 
shown in Fig. 1 compared with predictions from MFT 
[()]. As expected, the average activity levels of the simu- 
lated network coincide with mean-field predictions in the 
stable phases of low and high activity. However, in the 
transition region between the two phases, the activity 
levels from network simulations increase continuously in- 
stead of exhibiting the jumps and hysteresis predicted by 
mean-field theory. We note that the hysteresis in MFT [G] 
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Figure 2. (Color online) Firing polygraphs of the = 64 
network at 4 different values of synaptic strength w labeled 
in Fig. 1. For each value, only activities of 16 randomly chosen 
neurons are shown, with each occupying 1/16 of the height 
of the corresponding strip. The light and dark vertical line 
segments mark the firings due to noise and threshold crossing, 
respectively. The time on each strip of recording runs towards 
the left with a scale indicated by the one second time interval 
on each strip. Symbols (< and *) mark clusters of threshold 
firings as detailed in text. 



comes about when there are two stable states in the sys- 
tem: a quiet state with only noise-triggered firings that 
are insufRcient to ignite the entire system, and an active 
state in which firings are system- wide and self-sustaining. 
However, for a finite-size network having sufficient fluc- 
tuations and running for a sufficient time, the system 
can make transitions between the two states, leading to 
a single-valued mean activity level of the system show- 
ing no hysteresis. Comparing networks of different sizes, 
the rise of the average firing frequency in the transition 
region is steeper for larger networks since fluctuations in 
the total synaptic conductance decrease with an increase 
in the number of afferent synapses [6] rendering mean- 
field theory a better approximation. 

The firing pattern of the A'^ = 64 network is shown 
in Fig. 2 for 4 different values of synaptic input Kw as 
marked in Fig. 1. With computer simulations, we are 
able to distinguish firings due to threshold crossings of 
the membrane potentials (marked with dark vertical line 
segments) from noise triggered firings (marked with light 
vertical line segments). We can thus define a system to 
be active when there is any neuron with a membrane po- 
tential set to cross the firing threshold without the help 
of further noise events. Operationally, this is when the 
time-to-fire t\^^ calculated in step 1 of the event-driven al- 
gorithm on the preceding page is finite for any of the neu- 
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Figure 3. (Color online) Time-weighted duration distribution 
of sustaining threshold firing activities. The system is con- 
sidered active when any of the neurons is set to cross the 
firing threshold without further noise events. The duration is 
measured for each episode of the system staying continuously 
active. The oscillations in the distribution comes from the 
periodic nature of persistent firing activities. Values of the 
synaptic weight w (0.014, dashed line and 0.021, solid line) 
correspond to the locations B and C, respectively in Figs. 1 
and 2. 



rons. The switching between the quiet and active states 
of the network are evident on the strips B and C. This 
leads to clusters of threshold firings during which the sys- 
tem remains active. In the specific model we considered, 
there are two types of threshold firings clusters that can 
be identified from the polygraphs. The first represents 
bursting^ which are brief clusters (up to a few hundred 
milliseconds) of threshold firings of generally similar du- 
rations. These can be seen for all threshold activities on 
strip B and some on strip C as marked by the symbol 
in Fig. 2. The second type of threshold firing activ- 
ity is intermittently persisting activity, that can last for 
seconds, are of a wider range of durations and are seen 
in two time segments marked by the symbol "*" on strip 
C of Fig. 2. The two types of clusters of threshold fir- 
ings can also be identified from the duration distribution 
of self-sustained threshold firings shown in Fig. 3. The 
distinction between the two kinds of activities can be 
attributed to the short-term depression caused by the in- 
active state of synaptic transmitters in the TUM model: 
When threshold firings of the network are just ignited af- 
ter a quiet period, most transmitters are in the ready-to- 
rcleasc state and activities can propagate easily leading 
to a somewhat higher firing rate of the neurons in the 
beginning. However, as most transmitters arc deposited 
into the inactive state during repeated firings, the fir- 
ing frequency is reduced. At the low activity end of the 
transition region, the reduction of firing frequency typi- 
cally continues all the way to a cessation of any thresh- 
old firings, leading to the short bursting events of similar 
duration. Near the high activity end of the transition 
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Figure 4. (Color online) Scaling of "stopping time constant" 
estimated from the tail of active duration distribution in Fig. 3 
for a A*' = 32 static uniform network. The best estimate of 
uip = 0.0487 for the diverging point (center curve) is shown 
along with two slight deviations. The dashed line is to show 
the estimated slope of -2 for comparison. 

region, as the transmitters become inactive, the reduced 
firing frequency can settle to a stable value that can last 
for seconds before all threshold firings stop. From the 
simulation results such as what are shown in Fig. 3, the 
stopping of the threshold firings seem to follow Poisson 
statistics, that is, it can be characterized by a stopping 
rate r^^. The duration distribution of sustaining thresh- 
old firings decays exponentially for large duration with 
the time constant r that increases with w and appears 
to diverge at some w = Wp. For the TV = 32 network, we 
estimate Wp « 0.0487 and that the "stopping time con- 
stant" diverges roughly with the scaling r ~ {wp — w) 
as suggested in Fig. 4. We note that the scaling form 
of the stopping time remains the same when we simplify 
the TUM dynamics in our model by removing the inac- 
tive state (data not included here), suggesting at least 
some degree of universality. These results suggest that 
the divergence of the stopping time constant r signals 
that the system enters the phase of truly self-sustaining 
activity. 



V. DYNAMIC (PLASTIC) NETWORK 

For the current study, we focus on the stationary states 
of the system under the plasticity dynamics. A uniform 
network with synaptic weight Wj,i = w* for all synapses 
is used as the initial configuration, and simulations arc 
conducted at each data point for at least 2'^^ ms (« 25 
days) of simulated time to reach the stationary state. We 
concentrate our calculations on a fully-connected network 
of = 32 neurons with a high resolution in the plasticity 
parameter w* . For each w* , all measurements are aver- 
aged over an ensemble of 32 or more independent runs 
of different random sequences. The results of the calcu- 
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Figure 5. (Color online) Activity-plasticity(parameter) plot 
for a fully connected plastic network of A' = 32 neurons com- 
pared with the activity level of a static network of the same 
size with uniform weight w = w* . The plastic network is al- 
lowed to reach a stationary state by running for about 25 days 
of simulated time. All results are averaged over an ensemble 
of 32 or more runs of different random sequences. Insets of 
the plot show the normalized synaptic weight distributions on 
a log-log scale with a range of 2 (4) decades on the horizon- 
tal (vertical) axis at different points marked on the activity- 
plasticity curve. 

lations are summarized in Fig. 5, where we see that the 
firing frequency of a stationary plastic network coincides 
with that of a uniform static network [w = w* for all 
synapses) in the low-activity, noise-dominated regime as 
well as in the persistently active regime. In the transi- 
tion region between the noise-driven and self-sustaining 
regimes, the firing rate of the network is enhanced by the 
plasticity, and, while the behavior is not universal, with 
TUM synapses the amount of enhancement exhibits an 
interesting bimodal shape. 

In the insets of Fig. 5, we show the synaptic weight 
distribution, also averaged over the same ensemble, at 
each marked point along the activity-plasticity curve. In 
the two stable regimes, noise-dominated and persistently- 
active, the synaptic weights have a Gaussian distribution 
with a narrow width (within about 5% of w*) as predicted 
by MFT [()]. However, coincident with the bimodal en- 
hancement of firing frequency, the synaptic weight dis- 
tribution of the stationary plastic network shows dra- 
matic changes in the transition region along the activity- 
plasticity curve: Increasing the plasticity parameter w* 
from the noise-dominated low-activity regime, we see a 
discontinuous jump in the firing rate and sudden broad- 
ening of the weight distribution with tiered side peaks at 
w* w 0.017. The enhanced firing and broadened distribu- 
tion ease off after the jump, until these features become 
insignificant around = 0.04. At this point, with a 
slight increase of w*, two disconnected side peaks pop 
out in the synaptic weight distribution. However, such a 
splitting in the weight distribution is only accompanied 
by a gradual enhancement of average firing rate of the 
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plastic network. The two side peaks eventually merge 
back to the main peak with further increases of w* . This 
eventually returns the network to a uniform conformation 
with a narrow Gaussian synaptic-weight distribution be- 
fore the plasticity parameter reaches w* ~ 0.1. 

Instead of a simple increase of the Gaussian width, 
the observed broadening of the synaptic weight distribu- 
tion in the transition region just summarized (see Fig. 5) 
comes with structured side peaks, which signify emergent 
network conformations that we try to decipher and visu- 
alize with a pscudophysical approach in Section VI. For 
the = 32 and smaller networks, we are able to obtain a 
stationary state from the simulation run that is insensi- 
tive to the initial configuration at each data point (value 
oiw*). However, we do see runaway synaptic weights un- 
der the plasticity rules used for the iV = 48 and = 64 
networks in narrow, isolated intervals at w* w 0.0268 
and w* ~ 0.0194 respectively near the high-activity end 
of the transition region before the runs (2^^ms, ~25 days, 
of simulated time) are complete. For these instances, the 
diverging synaptic weights and firing frequencies of the 
driven neurons slow the simulation down to a near stand- 
still, and we can not always maintain a stationary system 
as we vary the plasticity parameter across the runaway 
point. We note that such runaways are commonly seen in 
models with Hebbian plasticity and arc typically "cured" 
with cutoffs in the range of synaptic weights. [13-16] We 
do not need such a device for the TV = 32 network results 
presented. For larger networks, when we do see such "iso- 
lated" runaways, we regard them as pathological [17, 18] 
and exclude them from the ensembles. 



VI. NETWORK "LAYOUT" 

Characterizing the structure of non-uniform networks 
is an active field of research [19-".^1]. Here we adopt a 
more intuitive and visual approach to extract any emer- 
gent structures of the stationary plastic network obtained 
above. Under our model setup, such structures reside in 
the framework of a fully-connected network of uniform 
connectivity, represented by the main peak in the synap- 
tic weight distribution. In this section, we attempt to 
uncover the topological meaning of the side peaks seen 
for the distribution in the transition region between the 
two stable phases of the system. 

In order to visualize the networks such that any emerg- 
ing features can be better revealed, we introduce a pscu- 
dophysical system of fictitious interacting particles in 
two dimensions representing the neurons, recognizing, of 
course, that the actual all-to-all network is not strictly 
two-dimensional. The nodes (pseudo particles) are imag- 
ined to be located on the plane, confined to the region 
near an arbitrary origin by an isotropic harmonic po- 
tential. The nodes interact with a (fictitious) pair-wise 
repulsive force that is determined by the synaptic weights 
between the neurons. To have the pseudo-particles rep- 
resenting neurons that are more strongly connected stay 



close to one another, we choose a repulsive pairwise force 

^^.J = 2 2 (10) 

to act on particle i due to j with e^^j = rj^i/\fj^i\ be- 
ing the unit vector along the fj^i = fi — Vj direction (in 
the two-dimensional pseudo-space). In fitting the entire 
layout of the pseudo-particles (representing neurons) to 
a fixed display area, the strength k of the confining har- 
monic potential is adjusted so that, in an equilibrium 
arrangement of the particles, the radial coordinate r of 
the particle farthest from the origin of the potential has 
magnitude equal to 1. Such a pscudophysical system 
of neurons is relaxed to reach a stable arrangement of 
neurons in which the net force on each pseudo particle 
(neuron) is zero. The relaxation is accomplished through 
a deterministic, over-damped dynamics in which the ve- 
locity of each node is proportional to the net force it 
experiences 

§.-...+1:4. (11) 

i 

The representative layouts for the marked data points 
in Fig. 5 are shown in Fig. 6, where we only show a frac- 
tion of the synapses with strongest weights along with the 
pseudo-positions to reduce obscurity. We note that the 
stable layouts reached are not unique for the system, and 
it is possible for the pseudo-particle system to be trapped 
in a mctastable configuration in which strongly connected 
neurons are hindered by other neurons and fail to get 
close to each other during the relaxation. We do not at- 
tempt to further relax the layout (by, e.g., perturbing 
the positions of the neurons) after a stable configuration 
is reached during the relaxation process. (We note that 
for larger systems where the "hindrance" can be more 
severe, one may consider a three or higher dimensional 
pscudophysical system to help alleviate the problem.) 

As revealed in layouts of Fig. 6, the jump in the firing 
rate near w* = 0.017 is accompanied by the formation 
of a path conformation in the network where synapses 
of stronger synaptic weights connect a sequence of neu- 
rons and facilitate the propagation of activities along the 
path. Such a path can form a closed loop near the jump 
but become less likely to do so as w* is increased. The 
path or loop structure is stable in the system such that 
once it forms, it will generally remain to the end of our 
simulation run at a fixed w* value. As we only present re- 
sults of simulations from uniform initial networks, there 
is no hysteresis loop in Fig. 5. However, the jump in the 
firing rate signals a first-order phase transition. While 
we have not attempted a systematic study of hysteresis 
effects with larger systems, we do see some degree of hys- 
teresis in smaller (A^ < 16) networks: If we use a loop 
conformation obtained in the stationary state of the net- 
work at a w* value above the jump as the initial condition 
for the simulation of the network at a w* value below the 
jump, the loop conformation can sometimes persist over 
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Figure 6. (Color online) Typical resultant layouts of neurons (represented by the dots) following the pseudophysical approach 
described in the text at the points marked on the activity-connectivity curve in Fig. 5. The number labeled below each layout 
is the corresponding value of the plasticity parameter w* . Only the strongest 31 synaptic connections (lines with arrows) are 
shown in each layout. 



the course of the entire simulation. With a path or loop 
conformation of the network, the tiered side peaks of the 
synaptic weight distribution in this area come from the 
skip-over synapses along the path, i.e., synapses that are 
along the direction of the path but bypass some of the 
neurons. 

Eventually, with the increase of w* , the path falls 
apart and the network returns to a uniform conforma- 
tion around w* ~ 0.04. The emergence of disjoint side 
peaks in the synaptic weight distribution with a slight 
increase of w* from 0.04 as seen in Fig. 5 is accompanied 
by the formation of a sink hub structure in the network 
where the high activity level of the hub neuron and its 
strong afferent synapses reinforce each other because of 
the positive-feedback nature of the synaptic plasticity. 
This type of hub structure is first seen at w* = 0.0411 in 
our simulations and can go in and out of existence several 
times during the course of our simulation run at a fixed 
w* . This partially explains why we only see a gradual 
increase of the average activity level when a hub start to 
form. The sink hub structure becomes more and more 
stable as w* increases until multiple hubs appear in the 
network around w* = 0.0414. The network returns to a 
uniform conformation through the increasing number of 
less and less significant hubs. Finally, at sufficiently large 
w* the network has reached a homogeneous, high-activity 
state. 

We can now try to understand the formation of "path" 
and "hub" conformations in different segments of the 
transition region separating homogeneous, stable low- 



and high-activity states of the network. Near the low- 
activity end of the transition region, the sequence of neu- 
rons along the path represents a firing order in a bout of 
threshold firings that is amplified and preserved by the 
STDP as the path of stronger synaptic weights. Espe- 
cially when such a path forms a loop, threshold firing 
activity can cycle through all neurons in the loop no mat- 
ter which one of them is triggered by the noise. However, 
near the high-activity regime, each neuron generally fires 
multiple times during an episode of sustaining threshold 
firings with a high frequency. The firing order of neurons 
loses its meaning because there are many ways of pair- 
ing up their spikes. In this area of the transition region, 
the main distinguishable feature of a neuron is simply 
its firing rate. In fact, the expectation value of a synap- 
tic weight w for a synapse in this range can be uniquely 
determined from the firing rates of its pre- and postsy- 
naptic neurons when the firing rates of the neurons are 
high enough (A > t^^,t^^) so that the actual timing 
of spikes becomes unimportant. Under such a condition, 
the stationary synaptic weight under the plasticity dy- 
namics (6) for, say, the j ^ i synapse, is given by 

w,,=w*^ (12) 

where Yi = (Yi) and A.; = {Si) are, respectively, the av- 
erage active transmitter fraction and mean firing rate of 
the neuron i. With Eq. (12), the condition 

Wt.jW-j^kWk,i = w*^ (13) 
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should hold for any neurons i ^ j ^ k in the network 
[22]. And, indeed, we have verified that the condition 
(13) is well satisfied by the hub conformations obtained 
in the simulations, but is apparently violated by the path 
conformations (data analysis not presented here). 

While it is tempting to correlate the layouts result- 
ing from our pseudophysical arrangement with configu- 
rations of physical neurons, we must caution that one 
can only view such a pseudophysical approach as an ar- 
bitrary and non-unique way of revealing network struc- 
tures. Furthermore, the structures found by the pseudo- 
physical approach do not represent any direct physical 
information contained in the parameters of the neuron 
network. Nonetheless, it helps to provide an intuitive 
"structural" understanding of the network conformation 
as we have shown above. 



VII. CONCLUSIONS 

The event-driven algorithm presented in Sec. Ill im- 
proves the accuracy and efficiency of the simulations al- 
lowing us to gain an intensive view of the phase space of 
a stationary plastic neural network. With a large number 
of parameters typical for a model of a biological system, 
we fix all but one, the plasticity parameter, with phys- 
iologically plausible values for our investigation. Also, 
to minimize conceptual complications, we study only the 
stationary properties of fully-connected networks driven 
by uniform, non-correlated Poisson noise. Our results 
show the network develops interesting structures in the 
transition region between the stable phases correspond- 
ing to low- and high-activity regimes. While our fully- 
connected network with a limited number of neurons is 
certainly inadequate to address the dynamics of a real 
brain, it can be a reasonable starting point for cultured 
networks consisting of hundreds of neurons with virtu- 
ally "all-to-all" interactions [23-27]. Our finding sug- 
gests that emergent structures in these networks are more 
likely to be seen when it is firing intermittently in the 
transition region perhaps during its development from 
a weakly-coupled system of neurons to a more strongly 
connected network. 

The characterization of network connectivity is an ac- 
tive field of research with established quantifications such 
as the "clustering coefficient" [28] and "modularity" [21]. 
However, in the current study, we adopt a simple visu- 
alization approach to gain a more intuitive view of the 
network structure. The resultant identification of the 
"path" and "hub" conformations explains or, more accu- 



rately, coincides, with the appearance of structured side 
peaks in the synaptic weight distribution and the ele- 
vated firing rate of the neurons. While these path and 
hub conformations are the simplest forms that can ap- 
pear in a connected network, the present work is the first 
to demonstrate that they can arise naturally in transi- 
tion regions of a neural network under STDP and driven 
only by noise. The connection- weight conformation (i.e., 
the distribution and correlations of synaptic weights) of a 
naturally occurring network is often strongly infiuenced 
by activity during its formation and maturation. How- 
ever, studies of network dynamics typically focus on be- 
havior of the network under given, static network topol- 
ogy or weight conformation. Advances in the understand- 
ing of the network plasticity are setting the stage for 
addressing how the observed conformation of these net- 
works can come about. Our study of a pure and simple 
network is just an initial step along this direction, and 
it leaves open questions of how variations of biological 
details that are fixed in our model and different network 
topology, geometry, or sparseness might affect the emer- 
gent structures and alter the dynamical behavior of the 
network. 

In addition, there remain pronounced and puzzling fea- 
tures observed in the current study that require further 
elucidation. For example, while both types of conforma- 
tions (hubs and paths) appear abruptly as the plasticity 
parameter w* is varied, the average firing rate increases 
with a discontinuous jump for the path conformation, 
but continuously for the hub. Also, the system does not 
"morph" from the path conformation to the hub confor- 
mation directly but, instead, returns to a uniform net- 
work before developing the hub structure. This feature 
might suggest a symmetry or "topological" difference be- 
tween the two types of conformations that prevents a di- 
rect transition from one conformation to the other. 

Finally, while we have only considered a stationary net- 
work driven by uncorrelatcd noise, arguably the ultimate 
"goal" of a plastic neural network is learning to serve spe- 
cific or varying functional requirements. On this regard, 
it will be of great interest to subject the system to more 
meaningful inputs and explore the change in the station- 
ary structures as well as the transient dynamics of the 
resultant network. 
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